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Matrices are said to behave as free non-commuting random variables if the action which governs 
their dynamics constrains only their eigenvalues, i.e. depends on traces of powers of individual ma- 
\^ • trices. The authors use recently developed mathematical techniques in combination with a standard 

0^ ' variational principle to formulate a new variational approach for matrix models. Approximate varia- 

0^ . tional solutions of interacting large-N matrix models are found using the free random matrices as the 

variational space. Several classes of classical and quantum mechanical matrix models with different 
^ types of interactions are considered and the variational solutions compared with exact Monte Carlo 

D ■ and analytical results. Impressive agreement is found in a majority of cases. 
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^ . I. INTRODUCTION 

The main motivation behind this work lay in developing approximate analytical tools for dealing with the long 
' range physics of Quantum Chromodynamics. The large-A^c approximation is an attractive candidate for this project. 
Initiated by 't Hooft, Ref. 0, it was immediately applied to solve large- A'^c 1-fl dimensional QCD, Refs. Q, 
\^ 1^ and led to extremely fruitful phcnomenological insights and models, Refs. [Q. The factorization property 

^\ , of correlators of gauge invariant operators suggested the idea of the master field - the gauge field configuration which 
dominates the path integral in the large- A^c limit, cf. Refs. Knowledge of the master field should allow to 

calculate any gauge invariant observable as if it were a classical non-fluctuating object. 
^ A concrete example of a master field was provided by the exact solution of a model which described the dynamics 

(classical as well as quantum mechanical) of a single hermitian N x N matrix, Ref. |9|. It was found that in the 
f~| : large- limit the Boltzmann integral (in the classical case) and the ground state properties (in the quantum case) 
^ of this model were completely determined by an ensemble of matrices with a frozen distribution of eigenvalues given 
. ^ , by a solution of a certain integral equation of the mean field type. The "angular" orientations of the matrices in 
• this ensemble, i.e. the unitary matrices W in the decomposition M — WmW^ with m — diag(mi, m2, toat) 
J-j ] were completely free, i.e. distributed with equal probability according to the unbiased U (N) group theoretical Haar 
measure. This work and, in its sequel, Ref. pO[ | gave rise to an intense interest in solutions of large- iV matrix models 
and made contact with the study of such models in other fields ofphysics, most notably nuclear physics, cf. Ref. |Tl| ] 
and, more recently, two-dimensional quantum gravity, cf. Ref. and quantum chaos, cf. Ref. fl^ . The general 
interest in large- matrix models has been behind a steady progress in understanding solutions of various versions 
and generalizations of such models as reported e.g. in Refs. Q, ||l5| , jl^, 

A new direction in this activity opened up when the physics community became acquainted with recent advances 
in the mathematical literature, associated with the work of Voiculescu, Ref. |p^ , on non-commutative probability 
theory, Ref. [ pO| . The statistical mechanics (classical or quantum mechanical) of matrix (non-commuting) variables 
or (gauge) fields is an example of probability theory of non-commutative variables. A concept which turned out to 
be especially useful for the study of large- Af matrix and gauge theories, cf. Refs. ||2l| and [Q, was the introduction 
in Ref. ||l9[| of the so-called free random variables. The simplest examples of such variables occur in independent 
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lavge-N matrix models, i.e. in the classical statistical mechanics of several iV x iV, N ^ oo hermitian matrices 
Mi, i — 1, . . . , D, with the probability distribution (the Boltzmann factor) given as const ■ exp{—N^ J^i ^i(-^i)) with 
Vi{Mi) = {1 / N)Tt Pi{Mi) and Pi polynomial functions. Typical "observables" in such a model are correlators 



D 



l(Tr (M,, . . . M,J) ^ / 1 Tr (M,, . . . M,J exp{X][i^. - N'V,{M,)]} f[ Df,[Mk] , (1) 

i=l fc=l 

with the integration measure 

N N 

Dfi[Mk] = n di^-^kh^ n d^<Mk)j. dIm(Mfe)7. , (2) 

7—1 7>i^— 1 

and Fi = — In J Z3/i[Mi] exp[— iV^yi(Afi)] . Were the Mi regular commuting variables, one would call them independent 
since {Mi-^ . . . Mi,,) = (Mi-^^) . . . (Mi^.) for any selection ii, . . . For the matrix model this is obviously not the case. 
The probability distribution depends only upon the eigenvalues rrii of the Mi . The expectation of products of the Mi 
on the other hand involves also the "angular" variables Wi in the decomposition Mi = WirriiW^ . These variables 
however are "free", meaning that they are weighted with the unbiased Haar U{N) measure. The integrals over the Wi 
can therefore be evaluated in a universal manner independent of what the potentials Vi are. The remaining integrals 
over the rrii can then be evaluated in the limit iV — > oo by a saddle point approximation as in Ref . [D . The remarkable 
fact is that this procedure can be formulated in terms of very general rules relating any correlator of independent 
matrices to a linear combination of products of various individual moments of these matrices, cf. [hol. These rules 



will be reviewed in Section [I B below. Any non-commutative variables distributed in such a manner that they satisfy 
these rules are termed "free random variables" . Refs. [ pl| and |22[ provide a review of the techniques which were 
developed to deal with such things as sums and products of free random variables, to define a certain operator algebra 
in a suitably defined Fock space, and further applications, cf. also |2j], pst . 

It has been possible to prove that a master field for any matrix model can be constructed on the basis of the 



operator algebra developed for free models, Ref. |22|. However, so far no dynamical principle has been advanced 
in order to carry out this construction explicitely for interacting matrix and field theories. This paper presents an 
attempt to make progress in this direction by using a variational method based on the inequality {e~^) > e~^"^^. This 
approach is fairly standard in statistical physics. The new element will be to use the free matrix models as the trial 
variational space. The details of the development are outlined in Section ^ below. Given a certain set of variables 
in terms of which the action is formulated, one can give a complete solution of the variational problem in the space 
of free matrix models by defining what will be called the free partner of the exact action of the model. Subsequently 
one observes however that the available variational freedom is even wider since one can first transform to a new set 
of variables. Mi — ^ Mi{Bj). In the new Bi it is again possible to define the free partner of the action and give the 
general variational solution. Thus the "art part" of the proposed variational method is the choice of the optimal Bi 
or in other words of the optimal parameters of the transformation Mi — > Mi[Bj). 

Unfortunately it is hard to give estimates of the accuracy of this variational approximation. As usual, one must 



gain experience from examples and comparison with exact solutions. In Section III, a number of such examples are 
considered starting with classical interacting matrices and then going to quantum models of first one and then two 
interacting quantum matrices. Only linear transformations Mi = J^j ^ijBj are considered and the variational results 
using free random variables are compared with exact Monte Carlo calculations and, where available, with analytic 
results. In almost all models considered, e.g. in the models of two interacting matrices with the action 

S = ^Tr [M^ + Ml + gMlMl) 

and various values of of many interacting matrices with 

D 



S = y ^Tr MfM'^ and S = — ^ V ^Tr MfM^ 



/4 



and in the quantum models of one and two interacting matrices, impressively good agreement is found with exact 
results. Some models are also presented for which the method (with the additional limitation to the linear M B 
transformations) failed completely. Presumably nonlinear transformations could improve the agreement but the 
authors have not investigated this. 
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II. VARIATIONAL METHOD FOR MATRIX MODELS 



In this paper, matrix models described by partition functions of the type 

D 

Z = e-^ ^ / exp(-iV25[Mi, ...,Md]) \[ (3) 

k=i 

are considered, with M^, k — 1, . . . ,Z?, representing hermitian N x N dimensional matrices, in the limit N ~f oo. 
A factor A'^^ has been pulled out in the exponent and in return, all matrix traces appearing in actions S will be 
accompanied by a factor Note that the term "action" is used very loosely here to denote essentially any weight 

appearing in the exponent of the Boltzmann factor in a partition function. Only actions S[Mi, . . . , Mjj] which are 
invariant under the global unitary "rotations" Mi UMiU^ will be considered. One can further distinguish between 
"classical" and "quantum" matrix models. By classical one means the models for which the above integral is just the 
ordinary integral with the integration measure given by (H) whereas in quantum models this is a path integral with 
Mi depending on (imaginary) time r, Milr), and the measure 

N N 

Af W Dii[Mk{T)]=N n X{d{Mu{T))^^ n dRe{Mu{T))^u dlm{Mk{T))^u ■ (4) 

0<T</3 0<r</37 = l 7>!^=1 

Here (3 is the inverse temperature and A/" is a normalization factor which in principle depends on At, but will not be 
important in the following. 

In both the classical and quantum cases one deals with integrals over many matrix variables and the real distinction 
will lie in the form of the action. In the classical case the action (the classical potential) is a real function usually 
represented as a sum of monomials in powers of the Mi, 

S.,^V{M,,...,Md)^^ J2 E C:^S^^::^n'^{Mj^^Mi:^...Mj^) , (5) 

with real coefficients C^^'^^^'^''"^^^. In the quantum case the action will be a sum of kinetic and potential terms 
1 ^ 1 

^Q^^ E E]^Tr[Af,(T + Ar)-M,(T)]V ^rV {M,{t), . . . , Md{t)) , (6) 

0<T</3 1=1 0<T<I3 

with V as in (||). 



A. The Variational Principle 

As was reviewed in the Introduction, only a limited number of matrix models are amenable to analytical solution 
so that various methods of approximation are called for. Here, the powers of the variational method, combined 
with the large-A'^ saddle point evaluation, will be explored. As is common with statistical integrals, the variational 
approximation will be developed on the basis of the inequality, cf . Ref. |26| , 

-{x) 



> e-^^> , (7) 

which expresses the convexity property of the exponential function - the average value of the exponential of a random 
variable x is greater than or equal to the exponential of the average provided x is real and the weights used in the 
averaging are positive. 

If one can find an action 5*0 [M] which is simple enough to allow the calculation of integrals like 
J n i?A*[M,]e-^^^"M and J D^,[M,]KiM,, . . . , A/,,)e-^^^°M 

k k 

for some simple fmictions K then one can approximate the integral as follows. One can represent 
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where (0)s means / Oexp{F - N'^S)Dfi{M). Using the inequahty @ one obtains 



\ I Sq 

or 

F<Fo + N^S-SQ)g^. (10) 

If Sq depends on adjustable parameters, one can find their optimal values by minimizing the right hand side of the 
above inequality. The resulting Sq can then be adopted as an approximation for S with which one should calculate 
the relevant correlators. As always with variational methods, it is hard to assess the accuracy of the approximation. 
Here, it will be tested in different classes of matrix models, using the so-called free matrix models (cf. below) as the 
variational space, by comparing (where available) with exact analytical results or with Monte Carlo simulations. 

B. Free Matrix Models 

Freeness is a concept recently developed in the mathematical literature associated with noncommuting random 
variables . While freeness is abstractly defined without reference to a specific realization of the random variables, 
in the context discussed here, namely in relation to large-A'^ hermitian multi-matrix models, it essentially means 
that the weight exp(i^ — N'^S) according to which the matrices are distributed contains no bias with respect to the 
relative SU{N oo) "orientations" of the matrices. In terms of the decomposition Mk — WkrukW^ with diagonal 
rrik and unitary Wk this means that the action S in the free models depends only on the rrifc. Consequently all the 
angular integrals in any given correlation function are weighted only with the group theoretical U{N) Haar measure 
and give universal results independent of what the action is. The remaining integrals over the uik in the large-TV 
limit are controlled by the saddle point of the effective action which includes the contribution from the Vandermonde 
determinants for each k, cf. examples below. As a consequence of these properties all correlation functions for free 
multi-matrix models are systematically accessible in terms of the moments o f individual matrices through the recently 
developed techniques, as will be elaborated upon below, cf. also Refs. Q, |2^, p^ . 

A family of noncommuting random variables Mi is called free if and only if 

(/i(M,J/2(Af,J . . . /„(Af,J) = whenever (/,(A^.J) = for aU j (11) 

where ij ^ ij+i (and in+i is identified with ii ). Note that the bracket notation introduced here following p9| ] 
means taking the expectation value of the trace of the expression inside, divided by N . In order to evaluate a general 
correlation function {pi{Mi^)p2{^i2) ■ ■ -PniMi^)), one can use 

((pi(Af,J - (pi(Af,J)) (P2(M^J - (P2(M,J)) . . . {pniM,J - (p„(Af,J))) = (12) 

and multiply out the terms in the expectation value. This recursively determines correlation functions of a certain 
order in terms of lower-order correlation functions until one arrives at an expression containing only the moments of 
the individual matrices. For two free matrices Afi and Af2 this algorithm yields, e.g., 

(piiMMA'h)) = (pM)) (p2(M2)) (13) 
(AfiAfaAfiAfa) = (Mf) (A'h)^ + (Af|) (Afi)^ - (Afi)^ (Afa)' (14) 

Building on this information, it has been possible, e.g., to determine, p9[ |, how the eigenvalue distributions of 
matrices which are free with respect to one another are additively and multiplicatively convoluted. 

Freeness in matrix models is a slightly more general concept than what is usually meant by the term "independent 
matrix model" , in which the potential governing the distribution of the matrices is a sum of terms for each individual 
matrix. Independent matrix models are certainly free; however, so is e.g. a two-matrix model with potential 

V{Mi , A//2) = ^Tr A/2 + i-Tr + ^Tr Tr Af| . (15) 

More generally, the potential of a free matrix model can be any function of traces over functions of the individual 
matrices^. Such more general free matrix models are almost as easy to solve as the independent ones. This is 



^To the authors' knowledge, this type of model was first considered in |27 
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best illustrated by an example. Consider the model described by (|T^). After the transformation of variables Mi = 
WiTTiiW^ , with diagonal nn and unitary Wi, the potential V{Mi,M2) depends only on the m^. Moreover, the 
transformation introduces the Vandermonde Jacobians which, when exponentiated, modify the potential into 



N 



/ nOO poo \ nc 

Y / p^^^"^ in\fi-fi'\p,{p')dfidfi'+ / p^{p)fi^) + n / 



Pi{n)n'^d^ 



In the large- A'^ limit, the distributions pi{p,) of the eigenvalues are "frozen" at a saddle point, so that the moments of 
the matrices have definite values 

^TrMf-xi lTrAf| = X2 (16) 

This means that the saddle point eigenvalue distributions Pi{^) can be solved for as if they were controlled by the 
potentials 

V{Mi)^{x2 + l)j^TrMl V{M2) = {xi + l)^TYMj . (17) 
For such quadratic potentials, one obtains the well-known semicircular distributions, cf. Ref. [[Ill, |9[ 



Pi(M) = ^J^-M^ P2(M) = ^J^-M^. (18) 
The constants xi and X2 can now be determined by the self-consistency conditions 



Xi 



X2 = {Mi) = / dpp^p2{p) 



^Xi^X2^ — - — . (19) 



Including such self-consistency conditions on selected moments of the matrices is the only additional step needed 
compared with solving independent matrix models. Knowledge of the eigenvalue distributions Pi{p) is sufficient to 
calculate any correlation function since the latter can be reduced using Eq. ( |l^ ) to a sum of products of individual 
moments, i.e. to the integrals (Mf ) = / p''pi{fi)dp. 

It was argued e.g. in Refs. |£l|, that a free behaviour of the matrix variables representing physical degrees of 
freedom may capture correctly the main features of the dynamics in the large- limit of nonabelian gauge theories. 
This hope and the wealth of available analytical techniques presented above render the space of free matrix models 
a good candidate for variational calculations. The main idea of the present approach is to choose the optimal action 
Sq in the variational principle (^) from among all actions which leave the matrices in terms of which S is given, free 
with respect to one another. 

C. General Variational Solution in the Space of Free Matrix Models 

In making a variational approximation based on free matrix models one must first decide which combinations of the 
original matrix variables Mi will be assumed to be free, i.e. in which variables to formulate the exact action S. One 
should realize that e.g. free Mi, M2 do not imply free Mi + M2, Mi — M2 nor vice versa (more about this below). 
The second step is to find the best trial action Sq which is free in the variables one has settled for. 

In this subsection, the solution of the second step is addressed. In other words, consider So to be free in terms of 
the original Mi. Under this limitation, it is possible to give a general recipe for finding the trial 5*0 which minimizes 
the right hand side of dTO). For this it is useful to introduce the following algorithm. Given a trace of any function of 
matrices Mi, P ^ (l/A^)Trp(A/i, . . . , Mu), write down its expectation value in terms of the moments of the individual 
matrices under the assumption that the matrices are free with respect to one another. This is easily accomplished for 



any multinomial expression using Eq. (12). Now rewrite the resulting expression with (M|) replaced by Tt M^ /N. 
This defines what will be called the "free partner" Pf of P (the "liberated" P, so to speak). For instance, the free 
partner of P = (l/7V)Tr (M1M2M1M2) is 

Pf = ^TrMf (TrM2)2 + ^^(TrMi)^ TrM| - ^(TrMi)^ (TrM2)2 , 
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cf. Eq. Evidently, the original function and its free partner have the same expectation values when evaluated 

in any free matrix model for the Mi (note that the expectation values of the products of traces appearing in the free 
partner factorize at large N). In general, the two expectation values will of course differ. Furthermore, if the free 
partner S'/ of any action S is used as a trial action, it will generate a free matrix model; hence the terminology. 

With the help of the above definitions, one can give the general solution of the minimization principle ( |To| ) in the 
space of matrix models which are free in terms of the variables Mi. The crucial observation is that, when 5*0 is to 
be chosen from among all free actions, the original action S and its free partner Sf lead to the same minimization 
problem, i.e. identical right hand sides of Eq. (10), since {S)so = {Sf)so by construction for any free So, as mentioned 
above. However, the minimization problem for Sf is trivial to solve: Since Sf itself describes a free matrix model, 
the best approximation to Sf in the space of free actions is Sq = Sf itself. Since the approximation to S is governed 
by the same minimization problem, Sq = Sf is at the same time the best free approximation to S. 

This result provides a general solution of the minimization problem in the space of matrix models which are free 
in terms of the original variables Mi. 

Note that the solution Sf oi the minimization problem defines a type of mean field approximation to the exact 
action S. Interaction terms in S which are sensitive to angular correlations between the different matrices are 
replaced by terms in which each matrix couples only to selected moments of the other matrices, i.e. to some mean 
properties independent of the angular orientations. This is entirely analogous to the standard development of the 
Hartree or Hartree-Fock mean field theories from a variational principle using the space of (properly symmetrized 
or antisymmetrized) product states. There, in second-quantized language, the combinations in which the mean field 
enters the single-particle Hamiltonian are determined by Wick's theorem; here, in complete analogy, the combinations 
in which the moments enter the free partner are determined by the free random variable axioms as encoded in Eq. 
(p^). In fact, (|l^) is nothing but Wick's theorem generalized from the usual bosonic and fermionic cases to the case 
of objects obeying the Cuntz algebra (this is the algebra obeyed by the Fock space operator representation of free 
variables, cf. Ref. |2l)). 



D. Variable Transformations 



The idea of the present variational approach was stated as choosing the optimal action 5*0 from among all matrix 
model actions which leave the matrices in terms of which S is given, free with respect to one another. In this there 
exists an additional freedom of choice which will now be exploited, namely the choice of matrix variables in terms of 
which S is given, i.e. in terms of which one develops the mean field approximation by constructing the free partner 
Sf. Up to now, only a fixed set of variables was considered - the original Mi. However, if one rewrites the action S 
in terms of other variables, 



S{{M,}) = Si{M,{{B,})}) = S{{B,}) 



(20) 



and again looks for the optimal free action, now in terms of the new variables Bj, approximating S, one may find 
a better approximation than the one found using the variables Mi. It should be emphasized that free Bj in general 
do not imply free Mi nor vice versa. Therefore, the accessible variational space using free matrix models is in fact 
much larger than was apparent in the last section. By allowing different sets of variables in which to formulate the 
problem to be approximated, one can even include into the variation models in which the original matrices Mi are not 
free with respect to one another. To formalize this idea, it is necessary to reexamine the derivation of the variational 
principle ([lO|). Carrying out a variable transformation Mi Mi{{Bj}), one has 



n^M[Bj]e-'^'^[^^l+'"^[^'l > exp \-{Fo + N^{S - {\nJ)/N^ - So 



/So, 



(21) 



I.e. 



F<Fo + N^S- (In J)/N^ - 5o)so 
with the Jacobian (which must be non-vanishing to have a legitimate change of variables) 

Dfi[M- 



J 



Dfi[B,] 



(22) 



(23) 



Note that this is the determinant of a DN^ x DN^ matrix, where D denotes the number of matrix variables. Now, 
one can use the theorem of the previous section with the immediate result that the optimal free choice of 5*0 in the 



6 



new variables is the free partner oi S — (In J)/7V^. Thus, while for any given set of variables, the optimal 5*0 is known, 
one may now try to further improve the approximation (i.e. diminish the right hand side of (^2|)) by trying different 
sets of variables leading to different S and J. 

In general, this is technically difhcult, since In J is not a multinomial expression, for which the algorithm of 
finding the free partner is easiest to carry out. Therefore, in this paper the treatment will be specialized to linear 
transformations , 

M, = ^c„S, (24) 
j 

Then J = (det(cy ))^ is independent of the matrix variables Bj. This greatly simplifies the treatment of (p2[); the 
Jacobian does not play any role in the determination of the optimal free action Sq for given new variables Bj and thus 
the best free approximation 5*0 to the given S is the free partner 5*^ of S. Using furthermore that {S — Sf)g^ — 0, 
one arrives at the residual minimization problem 

F < F/ - In J (25) 

where Ff is the free energy associated with Sf. It now remains to try different choices of Cy, leading to different S, 
Sf, and finally Ff, such as to minimize the right hand side of ([2^). 



E. Calculation of the Free Energy 

Before actually carrying out the minimization procedure, it is convenient to clarify a technical point concerning the 
actual evaluation of Ff given Sf. The best way to achieve this seems to be the following. By subtracting a constant 
Fref , independent of Cij , from both sides of (|25| ) , one arrives at the equivalent problem of minimizing Ff — \nJ — Fref ■ 
For Fref , one can choose a reference free energy, preferably of an independent matrix model with the same number of 
variables as S and whose potential Vref is a polynomial of the same degree as 5*. This specification is not necessary, 
but convenient from the point of view of keeping the expressions one deals with in practice as regular as possible. In 
many applications, VrefiBi) — {1/N)Tt ^^Bf is adequate. Then one can write 



Ff - Fref = - 1^^/ n D^i[B,]e-'''^f^^^ + ^''/ n ^^[^^ 

3 

= -Inl l[Df,[B, 



e 



-WV^^flB 



a=l 

- {aS f[B] + (l-a)Vr^j[B]) 

3 



a=0 







\a^\n[ \{D^l[B,]e-''"-^'f^''^°'^ (26) 



da 

3 

= N^ Cda{~Sf[B]-Vref[B])s. (27) 
Jo ' 

where 

S'f[B,a]=aSf[B] + {l~ a)Vref [B] (28) 

In practice, the matrix model described by S'j is only insignificantly harder to solve than the one described by Sf] 
on the other hand, expectation values like the ones contained in (|2^) are very easy to evaluate and thus constitute 
a convenient way of evaluating the right hand side of the variational principle (p5|). Thus, in applications, the best 
formulation of the variational problem seems to be to demand minimization of 

da {Sf[B] - Vref[B])s, - (In J)/7V2 (29) 

as a function of the transformation matrix dj, where, to recapitulate, S'f is given by ( [2^ ) and Sf is the free part- 
ner of S, which in turn arises from performing linear transformations on the original matrix variables Mi, i.e. 
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S{Mi) = S{cijBj) = S{Bj). Furthermore, J = (det(cij))''^^ , and Ke/[-B] is a conveniently chosen reference po- 
tential independent of Cy generating an independent matrix model. 

As a last remark, if one has succeeded in minimizing (^9|) by varying the transformation Cy , one will usually be 
interested in transforming back to the original variables; e.g., given the eigenvalue distributions of the matrices Bj, 
one may wish to extract the eigenvalue distributions of the matrices Mi. Here, the additive convolution techniques 
for free variables developed in ||l9|] find fruitful application, as will become clear in the examples to be treated further 
below. 



III. EXAMPLES 



A. Models Involving Two Classical Matrices 

The simplest type of interacting matrix model involves just two classical matrices; this provides a first testing 
ground for the variational approach developed in the previous sections. Without any deeper motivation, the actions 

Si = ^Tt [Alt + Mi + gMlMl) 

and 

^2 = ^Tr [Ml + M| + gMiM2MiM2) 

will be considered. The variational approximation will be compared mainly with Monte Carlo experiments, although 
there also already exist exact analytical solutions to some simple models of this type, see e.g. in the case of 82- 

The models considered in this section are invariant under the transformation Mi — > —Mi, and throughout the 
treatment it will be assumed for simplicity that, at the large- iV saddle point, {Mi) ~ 0. Of course, it constitutes no 
problem in principle to work without this assumption and to verify it from the solution; this would just unnecessarily 
complicate the notation. In more complicated models, it may become an interesting issue whether such a reflection 
symmetry can be broken spontaneously. 

Consider now the model described by Si . This model will be treated here in some detail to exhibit the new techniques 
in practice. Later examples will receive a more cursory treatment. Allowing for a general linear transformation of 
variables, 

M, = ^Q,S, (30) 


and inserting into to arrive at Si{Bj), one obtains for the free partner of the latter 

Sif - bi^TTBt + b2^TTBl + h2^TTBlTTBl (31) 

with the notations 

hi = cfi + 4i + gcjicji 

h = + C22 + gc?2C22 (32) 

612 = 4(CiiCi2 + C21C22) + 5(CnC22 + C?2C21 + 2C11C12C21C22) . 

Here the reflection symmetry (Bj) = has been assumed to carry over from (Mi) — 0. Choosing furthermore the 
reference action 

Vref = i:TT{Bt + Bl) (33) 



N' 



one has (cf. Eq. {^) 



S[f = (abi + 1 - a)j^TrBt + (^62 + 1 - a)^TrB| + abi2^TT Bf TrSf (34) 

Treating the last term in analogy to (E5|) in conjunction with ( p^ ) and (p^, one has to solve one-matrix problems 
for Bi and B2 with the moments (Bf) to be determined self-consistently. The explicit solution of the one-matrix 
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model with arbitrary symmetric quartic potential is listed in Appendix |^; the expressions for the second moments 
yield equations for the (Bf) which can be solved numerically for given a and Cy . Then {Sif — Vref) can be evaluated, 
again using the expressions for the moments given in Appendix Finally, evaluating the integral over a leads to 
the free energy, cf. Eq. (^9|), with the additional Jacobian piece (In J)/iV^ = ln(ciiC22 — C12C21). Minimizing the free 
energy numerically as a function of the Cij yields the following result, cf. also Fig. |^: For sufficiently small \g\, the 
original variables Mi and AI2 are the optimal ones for a free approximation, whereas for large g, it becomes favorable 
to switch variables by substituting Mi — Bi + B2 and M2 = i3i — -82- This change in behavior is of course simply 
induced by the competition between the different terms in the action. For sufficiently low g, the quartic part, which 
is already independent in the original variables Mi, M2, dominates. For larger g, it becomes favorable to accomodate 
the interaction term as much as possible. The competition is essentially linear in the sense that one optimizes the 
treatment of either one or the other part of the action and the switch between the optimal choices of variables Mi , 
M2 or Bi = {Ml + M2)/2, B2 — {Mi — M2)/2 happens suddenly at ga = 2; the optimal choice does not change 
continuously with g. 

The authors have considered the possibility that this transition in the variational approximation may signal a 
large-A'^ phase transition also in the exact solution of the model, visible e.g. as a discontinuity of the derivative of 
the exact free energy as a function of the coupling constant g. However, the Monte Carlo results do not show such a 
transition - the expectation of Tr M^Af| displays a smooth behaviour near gcr- Thus, the transition in the variational 
calculation seems to merely reflect a crossover between two physically different regimes rather than a genuine large- iV 
phase transition. A higher order phase transition can however not be ruled out on the basis of the data taken by the 
authors. 

Having ascertained the optimal choice of variables, one can proceed to give the corresponding approximations to 
the eigenvalue distributions. For 5 = ±1, the original variables Mi, M2 are best; the free partner of Si is 

Sif = ^Tr {Mf + Ml) ± ^Tr MjTr Af| (35) 

Using the formulae of Appendix ^ to solve these models, one obtains (assuming the eigenvalue distributions of Mi 
and M2 to be identical) the consistency conditions 

{Ml) = {Ml) =XM = ^ [txm{xIi + 18) + {xIj + 12)^/2] (36) 

solved by xm — 0.334 for g — 1 and xm — 0.476 for g = —1 (the relevant solution can be picked out by using positivity 
of Xm, etc.). The corresponding eigenvalue distributions are given by (cf. Appendix ^) 

Pm(A) = -{2X^ + XM + m^)y/w? - with = \{\l xl, + 12 -xm) (37) 

TT 3 V 

These distributions are plotted in Figs. ^ and ^, and compared with Monte Carlo results for 10 x 10 matrices. In 
the case oi g ~ 1, the dependence of the Monte Carlo results on the size of the matrices is exhibited in Fig. IJ. In 
this, as in all other cases tested, = 10 already seems to embody the large- A^ asymptotic bulk behavior rather well, 
up to the characteristic oscillations in the eigenvalue density induced by the eigenvalue repulsion. In this respect, no 
deviation was observed from the well-known finite- A^ phenomena observed in one-matrix models. 

For g = 4, on the other hand, one considers the best free approximation after the variable substitution Mi = B1+B2 
and M2 = Bi — B2, i.e. 

Sif = ^Tr {Bt + Bt) + ^Tr B^Tr B| (38) 



Here, the consistency condition becomes 



1 



{BD - {BD ^XB = ^ -Axb{IQxI + 27) + {lQx% + 18)^/^ (39) 



solved by xb — 0.131. Eqs. (|3^) and ( p^ ) imply that the Bi are both governed by the quartic potential (l/A^)Tr (6i?^ + 
8xBBf), and the corresponding eigenvalue distribution is (cf. Appendix ^) 

Pb{X) = -{12X^ + 8xB + 6m^)x/m^ - ^ with = 1-{J64:X% + 72- 8xb) (40) 

TT 18 » 

In order to obtain the eigenvalue distributions of the original variables Mi and M2 from this, one must use the 
free additive convolution techniques derived in In Appendix M, the convolution of two identical distributions 
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governed by an arbitrary symmetric quartic potential is discussed. The resulting eigenvalue distribution pMi = PM2 
is plotted in Fig. ^ along with the corresponding Monte Carlo results for N = 10. Also, the comparison with the free 
approximation in the original variables Afi, M2 is plotted. Evidently, allowing the choice of free variables to vary 
leads to a vastly improved approximation in this case. 

In exactly the same vein as the case <? = 4 one can treat the extreme case (corresponding, up to a rescaling, to 

Sl^d ^ . (41) 

Here, in the absence of any piece in the action which is free in the original variables and which even in the case of 
quite strong coupling (7 = 4 (see above) kept the behavior of the model reasonably regular, the disagreement between 
the variational calculation and the exact result is rather catastrophic, cf. Fig. ^. This is not hard to understand; in 
the case of the action S^^'^, the variables Mi, M2 can use configurations such as 

"'-(71) (-) 

where nii, TO2 are roughly N/2 x N/2 matrices, to preserve a low value of S^'^'^ while realizing very large eigenvalues 
(together with a concentration of very small ones). This is indeed what is seen in Fig. |[ Dominance of configurations 
such as (^ ) implies a strong angular correlation between the matrices; if one rotates the above matrices with respect 
to one another, e.g. such that the eigenvalues are arbitrarily reordered, S^^'^ will in general take a very large value. 
Such angular correlations are completely lost in the free partner 5'[|'^ = (l/iV^)TrMfTrM| and evidently even 
the best choice of linearly transformed variables Mi = Bi + B2 and M2 — Bi — B2, though incorporating some 
angular correlations between Mi and M2, cannot do much to alleviate this problem. It is however possible that a 
nonlinear change of variables exists for which the difficulty will be satisfactorily removed. The authors did not explore 
this direction. The dismal failure of the variational approximation in this case thus does not come as a complete 
surprise. Rather, it is gratifying to see how well the variational approximation still works in the previous, not quite 
as pathological, example described by 5*1 with the strong coupling of 5 = 4. One striking feature of the Monte Carlo 
eigenvalue distribution in the case of S{^'^ is the absence of the characteristic finite- oscillations observed in all 
other cases. This can be understood as follows: Usually, the potential confines the eigenvalues to a compact domain 
and, due to the eigenvalue repulsion, they tend to be equidistant. This (fiuctuating) lattice has a favored equilibrium 
position, leading to the oscillations in the eigenvalue density. However, when the action allows the eigenvalues to 
spread over the entire real axis, the correlations induced by the eigenvalue repulsion become insignificant and the 
oscillations in the eigenvalue density disappear^. 

Apart from plotting the eigenvalue distributions, which essentially contain the information about all the moments 
of the individual matrices, one can test mixed correlators for freeness properties. In particular, if the matrices Mi 
and M2 are free with respect to each other, then there are e.g. the following relations between correlators, following 
from the axioms of freeness: 

Ci = {Ml Mi) = {Ml) {Mi) = Cif 

C2 = {M^M^) = (M^) (Af4) = C2f (43) 
C3 = {MlMiMlMi) = {Mt){Mi)^ + {Ml)^{M^) - {MD^M^)^ = C3/ 

In the Monte Carlo calculations, both sides of these equations were sampled in order to give a measure for the deviation 
from freeness (i.e. strength of angular correlations) contained in the exact models. The results are tabulated in Table 1 
for different values of g m. Si. Note that for Si'^'^, the moments of the exact distribution diverge due to the evidently 
too slow fall-off of the eigenvalue distribution for large eigenvalues, making such a comparison impossible. 
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{M?) 


{Mf) 


Ci 


Cif 


C2 


C2/ 


C3 


^3/ 


-1 


.483 


.377 


.254 


.233 


.171 


.142 


.150 


.121 


1 


.336 


.197 


.106 


.113 


.0339 


.0387 


.0272 


.0317 


4 


.267 


.133 


.0585 


.0713 


.0116 


.0177 


.0085 


.0139 



Table 1 : Mixed correlators, Eqs. (p3|), for the model Si compared to the predictions for free Mj 



^The authors are indebted to U.-J. Wiese for this remark. 
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One observes how, for growing g, the interaction term introduces stronger angular correlations, evident in the 
stronger deviations from free predictions. Now, for sufficiently large g, e.g. g = A, the variational principle asserts 
that a better approximation is obtained with a new set of mutually free variables Bi. Using this new set of free 
variables, of course also the predictions for mixed correlators, derived in (^) for free Mi, change. This shall be 
illustrated here for the correlator (Af^M|). Inserting A/i — Bi + B2 and M2 = Bi — B2, using the axioms of freeness 



on Bi and B2 together with (Bi) = 0, and then substituting back Bi = (Mi 
obtains 



M2)/2 and B2 = (A/i - Af2)/2, one 



{M1M2M1M2) 
(M1M2) 



= {Bt + Bl) = l{Mf + Al 



2MIM2MIM2) 



1 



{Bt + Bl) - 4{Bl){Bl) - (M^M^) ~ -((Af^ + Mif 
{BD - {Bl) = 



4(AfiA//2)^) 



(44) 

(45) 
(46) 



where in the last line, identical distributions for the Bi were assumed. Putting these relations together yields {M1M2) 
in terms of the individual moments, 



{MlMi) = i(Af4 + Ml) - -^{Ml + Ml)^ 



(47) 



In the case g — A, the right hand side, using the individual moments from the Monte Carlo experiment, takes the 
value .0617, which indeed is closer to the exact value for (Af^ Ai^|) = Ci than the prediction using free variables Af^ 
quoted in Table 1. 

Consider now the model described by 

^2 = ^Tr (A/2 + A/2 + gA/iA/2A/iA/2) . 



After performing a general linear transformation as in 



one arrives at the free partner 



1 



S2f = {4,+4,)-TTB 
^22 



-■12 



4)^Tri?2. 



2 .2 1t^^4 



N' 



-TtB^ 



5C11C21 ^ 
1 



4gciiCi2C2iC22-^Tr B^Tr B2 



(48) 



Evidently, from the point of view of the variational approximation, this model is slightly simpler than the one 
considered above. Since the noninteracting part of S2 is invariant under linear transformations of the variables and 
subsequent construction of the free partner, up to trivial variable rescalings (note also that, as before, (Bi) = was 
used in deriving (^)), the choice of variables will always be such as to best accomodate the interacting term, for any 
coupling g. Carrying out the calculation of the free energy in complete analogy to the calculation for 5*1, this turns 
out to be the choice A/i — Bi + B2, M2 = i?i — i?2, corresponding to the free partner 



1 



N 

The corresponding consistency condition is 
(B!) = {BI)=xb ^ 



S2f = -rrTr {gBf + gB^ + 2Bf + 2BI) 



^TjBITtBI 



108g2 



(2 - AgxB)^ - 18g(2 - Agxe) + ((2 - Agxe)'' + 12.g) 



3/2 



(49) 



(50) 



From the resulting eigenvalue distributions of the matrices Bi one again obtains the eigenvalue distributions of the 
original variables Mi by additive convolution. The result for g = 2/5, leading to xb = .2522 in (^0|), is displayed in 
Fig. 1^, compared with the Monte Carlo result for TV = 40. Also, various correlators are tabulated in Table 2. 
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{Mf) 


{Mf) 


Ci 


Cif 


C2 


C2/ 


C3 




2/5 


.53 


.56 


.30 


.28 


.37 


.32 


.28 


.24 



Table 2 : Mixed correlators, Eqs. (k3), for the model S2, compared to the predictions for free Mi 
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The exact eigenvalue distribution turns out to be surprisingly similar to the semicircle of radius \/2 obtained at 
g = 0. This is well reproduced by the variational calculation despite a quite strong dependence on g of the individual 
coefficients in 5*2/. The effects of the g-dependence of the quadratic and quartic coefficients nearly cancel in the 
complete potential, leading to a quite stable eigenvalue distribution. 

Turning to the mixed correlators in Table 2, one does not observe a strong deviation from free behavior of the Mi. 
In this case, there is almost no room for improvement by assuming free Bi instead of Mi , leading to equation (|4^ ) for 
(MfM|). Indeed, one here obtains for the right hand side of ( ^ ) the value .28, the same as Ci/ in Table 2. 

The reader may wonder why only the relatively small value g = 2/5 was displayed in the comparison above. In 
fact, the model described by 5*2 becomes unstable for too large values of \g\. Note that the product of two hermitian 
matrices M1M2 is in general not hermitian; therefore, the interaction term, proportional to (MiAf2)^, is not bounded 
from below even for positive g. Strictly speaking, the model is unstable for all g, but for sufficiently small g, there 
is effectively a barrier posed by the independent part of the action preventing the system from spilling over to the 
region of large eigenvalues where the interaction will allow eigenvalues to grow without bound. A glimpse of this 
instability (for negative g, namely setting in at 5 = —4/9) was already given in Here, instead the instability 

at positive g was investigated in slightly more detail. According to Monte Carlo simulations at iV = 10, 20, and 40, 
the critical coupling lies in the interval g £ [0.4,0.45]. On the other hand, the variational approximation in the best 
basis, described by the action (|4S|), is stable up to g — 2. Beyond this point, there is no solution to the consistency 
condition ( |5C| ) . The authors also checked that there is no two-cut solution to (^9|) above 5 = 2 (note that (^^ is valid 
only for one-cut solutions). Thus, the variational approximation does qualitatively capture the phase structure of the 
model described by S2, albeit with a badly overestimated critical coupling. In this respect, the choice of variables 
Bi = (Ml -|-M2)/2 and B2 = {Mi — M2)/2 represents a drastic improvement over the original variables Mi. In terms 
of the latter, the free counterpart of S2 is simply ^2/ = (l/A^)Tr {M^ -|-M|), which entirely misses the unstable phase 
of the model (in fact, all the dependence on the coupling g). 

B. Classical Constant Matrices in Many Dimensions 

Consider a model of D classical matrices with an action of the type 

1 ^ 

SD = ^^Y.SiM.,M,) (51) 

Note the scaling of this action with D, which is necessary to retain a meaningful balance between it and the eigenvalue 
repulsion originating in the Haar measure of the Mi. Before proceeding, a comment is in order regarding the special 



form of (51) and the interpretation of D as the number of space-time dimensions. In (pl[), each matrix degree of 
freedom interacts with all others in the same way. If one wishes to attach physical meaning to the matrices as degrees 
of freedom living in space-time, this occurs most naturally when all the degrees of freedom are attached to the same 
space-time point. Such a model may become relevant if one has managed to decouple the different space-time points 
of a _D-dimensional large- field theory, e.g. as the result of an Eguchi-Kawai reduction (for a review, see Ref. p9|). 
Now, in view of the scaling of the individual terms in (|5^), one might hope that the detailed angular correlations 
between pairs of matrices become unimportant at large D and a free approximation becomes exact. This argument 
can be made rigorous e.g. in the Kazakov-Migdal model, where the leading large- -D behavior is described by a free 
matrix model |30[ . The proof however depends on the specific link variable structure of the Kazakov-Migdal model. 



In order to test this idea more generally, the models described by 

1 1^1 
■ -TiMfMf and = ^ _Tr M/M/ (52) 



5D2 



1 ^ 1 



were investigated in the present framework for different D. The reader is reminded that in the case of Sd2 for D — 2, 
the free variational approximation failed miserably (cf. previous section). In analyzing the models ( p2|) in the free 
variational approximation, the authors did not treat different choices of free variables; for simplicity, only the best 



^Presumably, the model described by Si will display a similar instability when the coupling g becomes too negative. This was 
not investigated further by the authors. 
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free approximation in the original variables Mi was considered. Then the best free approximation is easily derived: 
The free counterparts of (^2|) are 



S 



D2f 



D 



1 ^ 1 

— T — 



Tr M^Tr and 



D 



1 ^ 1 

— Y — 



Tr M/Tr M,^ 



and thus each variable Mi is described by the action 



Sd2{M,) = 2xD2^TTMf or SoiiM,) = 2.xi34^TrM-' 



(53) 



(54) 



independent of D, with xd2 ~ (Mi) and a;£)4 — {Mi) to be determined self-consistently. One obtains xd2 = 1/2 and 
XDi = l/\/8 and therefore the eigenvalue distributions 



P2(A) 



and p4(A) = - I V2\^ 

TT 




(55) 



Comparing this with the results of Monte Carlo experiments, cf. Figs. ^ and the correspondence is quite satisfying 
already at D = 10 (the simulations were carried out for TV = 10, some checks for iV = 20 did not reveal any significant 
differences). Also, comparing mixed correlators of two matrices (without loss of generality. Mi and M2) with free 
predictions (cf. equations (E3)) in the Monte Carlo experiments, one obtains good agreement for D = 10, cf. Table 3. 









Ci 


C'lf 


C2 


C2/ 


C3 


C3/ 


Sd2,D^3 


.62 


1.00 


.25 


.38 


.33 


1.00 


.18 


.62 


Sd2,D = 10 


.51 


.53 


.25 


.26 


.25 


.28 


.18 


.21 


Sd4, D = 3 


.49 


.46 


.20 


.25 


.12 


.21 


.09 


.16 


Sd4,D = 10 


.46 


.37 


.21 


.22 


.13 


.13 


.11 


.11 



Table 3 : Mixed correlators, Eq. (^3|), for the models (|5^), and predictions for free Mi 

As mentioned above, this agreement does not come unexpected. However, while in the cases considered here the 
Monte Carlo solution seemed to converge towards the free approximation at large D, it seems doubtful that this 
should be true in general. On the one hand, a solution with angular correlations, say with angles restricted to a 
certain sector as compared with some fixed reference matrix, costs a weight exp(— aiV^D) in the partition sum (there 
are 0{N^D) angles, and the coefficient a depends on the details of the angular distributions). On the other hand, 
the interaction itself is also scaled to be of order 0{N^D). Therefore, there should be a genuine competition between 
the two terms, the former favoring freeness and the latter preventing it. 



C. Quantum Mechanics of One Matrix 



For a wide class of potentials, the quantum mechanical ground state of a single large- iV matrix can be solved for 
exactly This is achieved by working directly with continuous (Euclidean) time and mapping the problem to a 
noninteracting fermion problem. If one discretizes time as in the original definition of the path integral, the problem 
can be interpreted in terms of many interacting classical matrices. In this form, the quantum mechanics of one matrix 
is hard to solve exactly; on the other hand, it becomes amenable to the variational method dealt with here. Interest 
in such a discretized form has arisen in recent applications to quantum gravity, pi] . For definiteness, the model 
described by 

Sqi = ^Tr (eV^«+i - M,)^ + gMt^ (56) 

will be considered here, where L is the number of time slices, and Mq is identified with M^, i.e. periodic boundary 
conditions are posited. Generically, there is a competition between two different interaction terms: The nearest- 
neighbor coupling and the potential term, which acts instantaneously. In order to best accomodate the former in 
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a free variational approximation, one will want to change matrix variables to a Fourier basis, which decouples the 
kinetic part of the action. On the other hand, the potential term will prefer the original variables. Thus, in the present 
treatment, not a general linear change of variables will be considered, but only the aforementioned two discrete choices; 
the variational minimization principle will decide which choice is preferred. The Fourier basis is defined as 



Mr 



(57) 



with [L is for simplicity taken to be even) 



U = 



1 ° 







1 


1 


1/V2 \ 


sinfcL/2-1 


sin ki 


l/x/2 


cos fci 


COSfcL/2-i 




sin2fci/2-i 


sin2fci 


l/x/2 


cos2fci 


COs2fci/2_i 


1/V2 


sin3fci,/2_i 


sin Sfci 


1/V2 


cosSfci 


COs3fcL/2_i 


-1/V2 


V sin(L - 1)A;l/2-i • 


• sin(L — l)fci 


1/^2 


cos(i — l)fci • 


• COs(i - l)fci/2_i 


-1/V2/ 



(58) 



Here, the wavevectors are kj = 27rj/L; note that U is orthogonal and therefore the Jacobian term in the free energy 
( p9| ) gives no contribution for the Fourier choice of variables. For convenience in notation, the index j labeling the 
variables Bj in (p7|), i.e. the columns in (|5^), will be taken to run from —L/2 + 1 to L/2. 

In the Fourier basis, the kinetic part of the action (Bq) is exactly diagonalized and coincides with its free partner 



j = -L/2+l 



L/2 

E 



On the other hand, the potential part in the new variables is 



^CjiiBi) — ^Tr ^ ^ UijBjUikBkUiiBiUimB,, 

i jklm 



(59) 



(60) 



Assuming as before (Bi) = 0, the only terms which will give a contribution to the free partner are ones in which the 
four indices j, fc, I, m all take the same value or pairwise two different values. Thus, one has the free partner 



j<k \ i / j \ i / 



(61) 



By explicit calculation, U^j is of order 0{1/L) and therefore the quartic term can be dropped for a large number 
of time slices L. On the other hand, again by explicit calculation. 



L-l 



E^^^^^ = 7 forallj,A: 



(62) 



i=0 



and therefore one finally has for the free partner of the action (^6|) in the Fourier variables B^, 

L/2 L/2 

SQifiB,)^2 il-co,h)-TrB^ + f ^ l^^^'TrB| 



i=-L/2+l 



i<j = -L/2+l 



(63) 



On the other hand, the free partner of (56) in the original variables is 

L-l 

N' 



1 

SqifiM.) = ^Tr ^ (2M,f + .gM,f ) 



(64) 



again assuming {Mi) = 0. Here, the potential term is treated exactly, whereas the nearest-neighbor coupling has been 
completely truncated. 
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The free models are now easily solved. For the purpose of calculating the free energies, the reference action 
VrefiCi) = (l/7V)Tr J2^ Cf ^ill bc uscd; then one considers (cf. Eq. @)) 



SqifiM^) - ^Tr ^ [(1 + a)M^ + agUf] 



(65) 



for the original variables and 



S'qiS^B,) = ^(1 + « - 2acosfc.)^Tri?f + ^Y. j^TrS.^TrS^ 

i<3 



N' 

for the Fourier variables. In the latter, introducing the abbreviation 



(66) 



(67) 



one has semicircular distributions for the matrix variables Bi of radius squared 

2 



* 1 + a — 2q; cos ki + agx b 
implying second moments {Bf) = r|/4, which self-consistently determines xb'- 



xb 



1 



+ a — 2a 008(2^1 / L) + agxB 



L — *oo 



1/2 



dk- 



1^2 1 + a + agx b — 2a cos 27rfc 
1 



y/{a + agxB + 1)^ — 4a^ 



(68) 



(69) 



The free energy is then 



L — *oo 



x^{a + agxB + 1)^ — 4a^a;^ — 1 = 



, 1 2 - 2cos(27ri/L) + .ga:B - 1 
da y ; 

^ 2 1 + a - 2q; cos(27ri /L) + agxB 



da 



„ 1 2- 2cos27rfc + .ga:i3 - 1 
dk —- 



= L da — 
2a 



1/2 2 1 + a — 2a COS 2nk + agx B 
1 



1 - 



\/{a + agxB + 1)^ — 4a^ 



(70) 



(71) 



which must be calculated numerically (remember that at every a, is determined by the equation (|70|)). 

On the other handjin the decoupled basis, one has the quartic action (^5|), identical and independent for all the 
Mi. Using Appendix |A|, the resulting moments are 

^ ^ '{1 + a)^ - 18ag(l + a) + ((1 + af + Uagf^ 



108a2.g2 
1 



216aV 

and, inserting in the free energy. 



(1 + af + 18ag(l + af + 54aV - (1 + ")((! + af + I2ag) 



3/2 



(72) 
(73) 



Fm = / da2L{M^) + Lg(M^) - L{M^) 



108g2 



-4-36g 



8lg' 



+ (4 + 305)v/l + 3.g + 54.g2 ln(l + ^1 + 3.g) 



(74) 
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Numerically, Fm = Fb at g = 1.2306. For smaller g, the Fourier basis provides the better approximation, for larger 
g the original choice of variables does. It has been argued in that matrix models such as the one considered 
here exhibit a Kosterlitz-Thouless transition at some value of the lattice spacing separating the two regimes where 
the matrix chain behaves more like a collection of independent sites and the one where it supports spin wave type 
of modes. The variational approach presented here, while giving no details of the transition region, is capable of 
capturing these two different possible regimes. 

It remains to give the eigenvalue distributions in the two regimes. In the original variables, described by the action 
SQif{Mi), cf. equation (|64|), one immediately has, using Appendix 



Pm{X) - -(2.gA2 + 2 + m^g)Vr 



A2 



with m 



3.9 



(v/TT3g-l) 



(75) 



for all the variables M,. In the Fourier case, one must convolute the semicircular distributions of the Bi to obtain the 
distributions of the Mi. It is well known JlOt that the semicircular distributions are closed under additive convolution 
and that the square radii add up to the square radius of the resulting semicircle. I.e., 



Mi 



2 2 



(76) 



Inserting (p8[) for a — 1 and (pq), one obtains identical radii for all the Mi 



L/2 



''"Mi 



^ 1=^/2+1 ^ " ^ cos(2^j7L) + gxB 



1/2 



1/2 2 — 2 cos2TTk + gxB 
2 



(77) 



\/^gxB'+'g^x^ 

2XB 



where in the last hne it has been used that xb = (2/i)(^- Bf) solves (cf. ( [70|) for a = 1) 

Agxl + g^xj, = 1 
The corresponding semicircular distribution for the Mi is then 

Pm{X) = —^J2XB~\^ 

XBT^ 



(78) 



(79) 



The variational approximation derived above should now be compared to Monte Carlo experiments and the exact 
results known about this model. In Figs, p^p^ Monte Carlo results for L = 10 and = 10 are exhibited together 
with variational results at different couplings g. Checks with iV = 20 and i = 20 revealed no significant differences. 
One indeed observes, as already borne out by the calculation of the free energies, that the Fourier basis is preferable 
in a weak potential and the original basis is favored by a strong potential. This is corroborated by the Monte Carlo 
evaluation of correlators between matrices at neighboring time slices, cf. Table 4. 
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Ci 


Cif 


C2 


C2/ 






.1 


.55 


.57 


.41 


.31 


.59 


.33 


.52 


.26 


1 


.24 


.11 


.066 


.059 


.015 


.012 


.014 


.0093 


10 


.098 


.017 


.0099 


.0097 


.00031 


.00029 


.00025 


.00024 



Table 4 : Mixed correlators, Eqs. (^3|), between matrices Mi at neighboring time slices and predictions for free Mi 

At large coupling, neighboring matrices are to a very good approximation free, meaning the original decoupled 
basis gives a good description. By contrast, at low g, there are significant angular correlations induced by the nearest- 
neighbor coupling. These angular correlations are washed out by the fluctuations when one considers matrices far 
apart on the temporal lattice, cf. Table 5. There is virtually no deviation from free behavior between distant matrices 
even at low values of g. 
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(Mi) 






Cif 


C2 


C2/ 


C3 




.1 


.55 


.57 


.31 


.31 


.35 


.33 


.28 


.26 


1 


.24 


.11 


.059 


.059 


.012 


.012 


.0093 


.0093 


10 


.098 


.017 


.0097 


.0097 


.00029 


.00029 


.00024 


.00024 



Table 5 : Mixed correlators, Eqs. ([13[), between matrices Mi separated by half of the total length of the temporal 

lattice, and predictions for free Afj 

In order to make contact with physical (still Euclidean) time, one should write the discretized action of the quantum 
one-matrix model as 



1 



(M' 1 - M[f L 



(80) 



where T is now the length of the circle in time direction. Rescaling Af/ — Mi y/2T/L, one regains the form (^6|) of 
the action, with the identification g — AgT^/L^. Therefore, g behaves like the third power of the lattice spacing. 

In the continuum limit, the ground state distribution of the eigenvalues of the quantum mechanical one-matrix 
model has been found analytically in This can now be compared to the variational solution. According to the 
correspondence established above, the limit of vanishing lattice spacing corresponds to vanishing g; therefore the 
Fourier basis is the preferable one. The eigenvalue distribution for the rescaled variables Af/ is (cf. equation dTE^ 



1 



2T xbtt 



(81) 



and is determined by (cf. equation (|7q)) 



which for large L is solved by 



Therefore 



2^3 



16g JqXb 



xb 



0(1) 



(82) 



(83) 



PM' (A) = - - 




(84) 



Note that there is no T-dependence left in the continuum limit. By contrast, the exact solution is 

8 



1 



p^?r*(A) = -V2e-2gA4 with 



SlTT^ff 



1 



r(i/4) 



(85) 



determined by the normalization condition. The two solutions are compared for 5 = 1 in Fig. 



D. Quantum Mechanics of Two Matrices 

While the case of one quantum mechanical matrix discussed in the previous section is amenable to exact treatment 
due to its special form, an analogous solution is not possible once more than one quantum mechanical matrix is 
involved. On the other hand, the approximate variational approach, using discretized time, is easily generalized to the 
case of more than one matrix. Here, a commutator-type interaction will be considered, as it is especially interesting 
from the point of view of applications to Yang-Mills theories: 

Sq2 = ^Tr |^^(Afi,,+i - M^,,f + {M2.,+i - M^.^f + g{i[Ml,^, M2Af^ (86) 
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where the first index labels the two different matrix variables and the second one the different time slices. Again, only 
the original choice of variables will be compared to a Fourier basis (cf. (|5^)) using the variational criterion. 

Note that the action Sq2 is invariant under simultaneous shifts of all the matrices by a multiple of the unit matrix. 
This trivial freedom should be removed in order to obtain stable Monte Carlo results; otherwise, the system just 
performs a random walk in the trace of the matrices. Thus, the Mn^i here are constrained to be traceless, i.e. the 
condition (Mn^i) — 0, assumed to be realized dynamically in prior examples, is enforced by hand. 

The free partner to (|8^) in the original basis is 



SQ2f (M„,,) = E ]^Tr {2Ml + 2MI,) + ^Tr Af^ ,Tr M|,, (87) 



i=0 



(using (Mi_i) — (M2,i) — 0). On the other hand, in the Fourier basis 

M„,, = Ec/yB„,, (88) 

3 

with U as in (p8|), the kinetic part becomes 

2 ^'"^ 

~S^42iBn,i)^~S^ci2f{Bn,) = J^Tr (1 - COS fc, ) (B?,, + ) (89) 

j = -L/2+l 

The potential part on the other hand becomes 

- ^Tr ^[Mi,,, M2,.]2 = ^Tr ^ ^ U,jU,kUuU,^{Bi^jBi^kB2,iB2,m - Bi^,B2,kBi^iB2,m) (90) 

i i jklm 

Using again the assumption (Bn.i) = 0, the second term in the round brackets on the right hand side never gives a 
contribution to the free partner, whereas the first term only gives a contribution if j ~ k and / = m. Thus one has 

S^Q2f{Bn,^) = ^ ^ J] C/., C/y C/.fc C/,fc = ^ E ^r (91) 

j^k i j,k 

Now one can again easily solve the free models. For the purpose of calculating the free energies, the reference action 

VrefiCn,^) = ^Tr ^(^1',. + ^^2,.) (92) 

i 

is convenient; then one has to solve the models 

^Q2/(M„,.) = J2 ^r^' (^^'^ + + ^Tr Mi^.^TrM^, (93) 

i=0 

for the original variables M„ ^ and 

SQ2f{Bn,^) ^ J^^l + a - 2acosfc,)^Tr {Bl^ + B^ + ^ E TrS^^.TrSf,, (94) 

j j,k 

for the Fourier variables Bn,i. Starting with (|9^), abbreviating {M'^^) = xm, one obtains identical semicircular 
distributions for all matrix variables with radius squared 

(95) 



a + 2agxM + 1 

whereupon the consistency condition determining xm becomes 



^^^ = T^L + 2a5XM + l ^^^^ 
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solved by 

XM = -^(V(l + a)^ + 4ga- (1 + a)) 

The free energy then is given by 

Fm ^ L I da (4a; Af + Sgx^j - 2xm) 
Jo 

On the other hand, in the Fourier case (p^, abbreviating 



(97) 



(98) 



(99) 



one notices that the variables Bi j and are controlled by exactly the same potential as the variables Bj in the case 
of the quantum mechanical one-matrix model discussed in the previous section. In particular, one has the consistency 
condition 




and the second moments 



The free energy now is 



x'^{a + agxB + 1)^ — Aa^x^g — 1 = 



(B^ ) = (B^ ) = - - 

^ ^ 2 1 + a - 2a cos kj + agxs 



L da 
'0 



2xB 



L/2 

E r 

j = -L/2+l 



2 cos{2tt j/L) 



9x% 



a — 2a cos{2ii j / L) + agxs 



xb 



(100) 
(101) 

(102) 



Numerical evaluation of the two free energies yields Fb < Fm for all i.e. the Fourier basis is always preferable over 
the original basis. While this is to be expected at small g, it is not necessarily an indication that the Fourier basis 
represents a particularly good approximation for large g as well: Rather, working in the original basis still truncates 
the interaction term rather badly at every time slice (cf. equation (^)). This is different from the one-matrix case, 
where assuming the original variables to be free merely implied decoupling variables at different time slices, which 
certainly should become exact for large coupling g. In the two-matrix case, by contrast, there is still an additional 
truncation at every time slice which prevents convergence to the exact problem at large g. In light of this, it is not 
so surprising that the Fourier basis is preferable for all g. Nevertheless, in Figs. |l4-16, this basis appears to provide 
quite a satisfactory description even for g — 10; the agreement with Monte Carlo results {N — 10 and L — 10 were 
used in the latter) is not much worse at g = 10 than at g ^ 1/10. 

Also in the correlators, cf. Tables 6 and 7, one observes that the angular correlations involving the same matrix 
variable at different time slices are stronger than the angular correlations between the two different matrix variables 
at a fixed time, up to quite high g. At g — 10, these two correlations seem to be roughly equally strong. This again 
corroborates the result that using the Fourier basis, which treats the nearest-neighbor coupling exactly, provides a 
good approximation for a large range of g. 



9 


{Mid 


{Mid 




Gif 


C2 


C2/ 


C3 


Czf 


.1 


.86 


1.5 


1.2 


.75 


5.6 


2.3 


5.1 


1.7 


1 


.38 


.29 


.19 


.14 


.16 


.084 


.13 


.063 


10 


.18 


.068 


.038 


.033 


.0064 


.0047 


.005 


.0034 



Table 6 : Mixed correlators, Eqs. 



j), between matrices Af„^i for one fixed n at neighboring time slices i and 
predictions for free Mn.i 
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{M?J 


{Ml,) 




Cif 


C2 


C2/ 


C3 


C3/ 


.1 


.86 


1.5 


.66 


.75 


1.7 


2.3 


1.3 


1.7 


1 


.38 


.29 


.12 


.14 


.059 


.084 


.043 


.063 


10 


.18 


.068 


.027 


.033 


.0027 


.0047 


.002 


.0034 
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Table 7 : Mixed correlators between the matrices Mi i and M2 i at a fixed time slice i and predictions for free Mi i, 

M2,.. 



Ultimately, though, one would expect a local treatment, which assumes different time slices to be free with respect 
to one another, to provide better agreement at high g if it manages to well approximate the commutator interaction. 
This is precisely what is not achieved if the original variables are assumed to be free. On the contrary, a nontrivial 
mixing of the variables Mi^i and M2,i at every time slice is necessary such as to better accomodate the commutator- 
type interaction in (^6|). Note that linear combinations of Mi^i and M2,i will not effect any improvement, since the 
commutator is invariant under such linear transformations (up to trivial rescalings of the variables). It is disappoint- 
ing that precisely this phenomenologically very interesting type of interaction term seems particularly resistent to 
free approximation combined with linear transformations. Here, it seems that using (technically more complicated) 
nonlinear variable transformations must be contemplated in order to capture the essential angular correlations. Note 
also that the commutator term in the action on its own does not confine the eigenvalue distributions of the involved 
variables to a compact support, even if one enforces tracelessness; the eigenvalues can become arbitrarily large if the 
matrices are located in the regions of configuration space where they commute. This pathology however seems to dis- 
appear (according to Monte Carlo experiments carried out by the authors) when more than two matrices are involved 
(i.e. a higher-dimensional Yang-Mills type of action). Presumably there, the regions of configuration space where all 
the matrices commute are too small compared with the whole space such that entropy suppresses configurations with 
arbitrarily large eigenvalues. 



IV. SUMMARY 



In this work, a new variational approach to interacting large- multi-matrix models was developed. The partition 
function was approximated using the variational principle F < Fq + {S — So) So with So initially taken from the space 
of all matrix models which are free in the set of variables the original interacting problem is given in. It turned out 
to be possible to give a general solution to this variational problem for a fixed set of matrix variables in terms of the 
concept of the "free partner" introduced by the authors. The free partner defines a type of mean field approximation 
to the original action and is constructed using the axioms of freeness in a fashion analogous to the construction of 
more conventional mean field theories for fermionic or bosonic particles using Wick's theorem. The freeness axioms 
constitute the analog of Wick's theorem for objects obeying the Cuntz algebra, which is the algebra obeyed by the 
Fock space operator representation of free random variables pl| . 

However, the variational approach presented here goes beyond simply acting as a device to derive a mean field 
theory. Above, the variational space was characterized as being the space of all matrix models which are free in 
the set of variables the original interacting problem is given in. This implies that one can considerably enlarge the 
variational space by first allowing for a change of variables in the original problem and only then varying over all 
matrix models which are free in the variables one has settled for. In this way one includes into the variation models 
which are not free in the original variables. Every set of variables defines a different mean field theory given by the 
corresponding free partner. The variational principle not only allows to derive these mean field theories, but further 
allows to decide which out of a set of mean field theories provides the best approximation to the exact problem. 

A number of examples were considered in order to assess the accuracy of the variational method; the variational 
results were compared with Monte Carlo simulations as well as analytical results, where available. To begin with, 
classical two-matrix models described by actions with at most quartic terms were investigated allowing for a general 
linear transformation of matrix variables. Impressive agreement was observed, except in the case of the action S^'^'^, Eq. 
(^). The authors did not attempt to improve the approximation by allowing for nonlinear variable transformations, 
which may relieve the discrepancy. Classical models with a large number of matrices were considered in which all 
matrices pairwise interact in the same way. Such models e.g. become relevant when one has managed to decouple the 
different space-time points of a large- A'^ matrix field theory in a large number of dimensions D, e.g. in the framework 
of an Eguchi-Kawai reduction. Among the models considered was one which reduces to for D ^ 2. It turned out 
that the quality of the variational approximation using only the original set of variables improves as D rises. 

The variational approximation was also considered for the problems of a single and two coupled matrix chains 
simulating path integrals for one and two quantum matrices. Two discrete choices of free variables, namely the 
original variables and a Fourier basis, were considered. In the single chain case, as expected, for small (large) lattice 
spacing, the Fourier (original) basis provided the better approximation. In the continuum limit of the one-matrix case, 
good agreement was achieved with the available analytical solution. For two matrix chains coupled via a commutator- 
type interaction, the Fourier basis provided a better approximation for all couplings; it also compares quite favorably 
with the exact Monte Carlo results for a large range of coupling constants. However, for very large coupling, one would 
expect a local basis, constructed such as to capture the angular correlations introduced by the commutator interaction. 
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to ultimately be more appropriate. Here, again it seemed that technically more involved nonlinear transformations 
must be contemplated in order to achieve further progress. 



APPENDIX A: SOLUTION OF THE ONE-MATRIX MODEL IN A SYMMETRIC QUARTIC 

POTENTIAL 

Following [|, H] the one-cut eigenvalue distribution p resulting from a hermitian large- iV one-matrix model in a 
quartic potential V{x) = 042;'* -I- a2X^ is most easily found by considering the resolvent, 



G{z) 



— 00 



z- A 



(Al) 



where Q{z) is a polynomial. Q{z), along with m^, which defines the support of p, is fixed by comparing coefficients 
in the asymptotic behavior of G, 



n—1 

where c„ is the n-th moment of the distribution p. One obtains 

1 



804 



- — i\ al + 12a4 - 02) 



Q{z) — 2aiZ^ -1-02-1- m^a4 rr, 
Then one can also immediately read off 

p{\) = -- lmG{X + ie) = -Q{\)^/m^ - 



(where e — > 0), and by expanding G in 1/z, 

1 



C2 



C4 



108al 
1 



180402 -f (oj -|- I204) 



3/2 



216a| 

When 04 = 0, these expressions simplify to 



02 -l- 180304 + 54o| — 02(02 + 1204) 



3/2 



P(A) 

C2 



1 



a2VW 



02) 



A2 



C4 = 



1 

20^ 
1 

2^ 



i.e. the second moment of the semicircular distribution is one fourth of its square radius. 

In complete analogy, one can derive the two-cut solution in the quartic potential. Then, G has the form 



G{z) = 2o4Z'^ -I- 022: — qyj (z^ — m\)(z'^ — m?_)z'^ 
from which one again obtains by comparing coefficients in the asymptotic expansion 

2 0.2 



q — 2o4 



± 



1 



204 ^04 

One can immediately read off the eigenvalue distribution, and the second moment e.g. becomes 

O4 y a 42 24 6\ 

C2 — — (m^ — m_|_m_ — m_^_ni_ + m_) 



(A2) 

(A3) 

(A4) 

(A5) 
(A6) 

(A7) 
(A8) 
(A9) 



(AlO) 



(All) 



(A12) 
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APPENDIX B: FREE ADDITIVE CONVOLUTION OF TWO IDENTICAL DISTRIBUTIONS 

GOVERNED BY A QUARTIC POTENTIAL 



In order to additively convolute the eigenvalue distributions of two matrices which are free with respect to one 
another, i.e. to obtain given and psj, one should go to the corresponding R-transforms These are 

defined as follows: Find the resolvents Gbi and as defined in Appendix 0; then the R-transforms are given by 
R{y) = G~^{y) — 1/y (here denotes the inverse function, G~^{G{x)) = x). The R-transforms behave additively, 
i.e. Rbi+B2 — R-Bi + Rb2- this sense they are the free analogs of the logarithm of the Fourier transform for 
probability distributions of ordinary commuting random variables. 

In the case of matrices governed by a quartic potential V{x) = a^x'^ + a-ix^ ^ the resolvent has already been given 
explicitly in Appendix]^; denoting G(z) = y, one has that z = G~^(jj) satisfies 

y = 2a4Z^ + a2Z — (2a4Z^ + a2 + rr?a^\/ 2^ — (Bl) 

where iv? was defined in Appendix Isolating the square root term and squaring the resulting equation leads to a 
third order equation for z. 



4a4yz — Aa^z + 2a2yz — y 



m^(a2 + rn^aiY' ~ 



(B2) 



It is not necessary to solve this explicitly, since the case of interest here is the case of two identical distributions, i.e. 
identical R-transforms. Then one has 



W : 



GbI+b^^v) 



RB1+B2 



RB,+RB2 + -^G-s\{y) + G-sl{y) 



implying 



2z- 



(B3) 



(B4) 



However, the equation satisfied by z is known, equation (B2), which by inserting (84) leads to the equation satisfied 
by w = Gg^^g^{y). Since one is actually interested in y = Gbi+B2{'^)^ can solve directly for y by rearranging 
(B2) in conjunction with (B4) as a polynomial in y, 



2y^ — y^{a4W^ -f 2a2w) — j/^(a4W^ -I- 2a2 — 2m^(a2 + m^a4)^) -|- ya^w + a4 = 



(B5) 



Rather than obtaining a cumbersome analytical solution of (B5), the relevant solution was tracked numerically. 
This is best done in small steps in along the positive real w-axis starting from the known asymptotic behavior 
y = \/w -\- 0{w~^) and using the positivity of the eigenvalue distribution psi+s^- The latter is essentially given by 
the imaginary part of y, namely PB1+B2W — ~ImGBi+B2(A -|- ze)/7r. 
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FIG. 1. Free energies of free variational approximations to Si as a function of the couphng constant g for different sets of 
linearly transformed variables Mi = ^ . CijBj, where Cij is orthogonal, Cij = afj cos cp + ajj sin cj> (and a'^ are the Pauli matrices). 
Note how all trajectories cross at Qa — 2; (j> = provides the best approximation below ga and (j> — 7r/4 above ga- Trajectories 
generated by nonorthogonal Cij always lie higher and are not displayed. 

FIG. 2. Eigenvalue distribution p(A) generated by Si a.t g = —1. Data points; Monte Carlo result for A'^ = 10. Solid line: 
Best free variational approximation. 

FIG. 3. Eigenvalue distribution p(A) generated by Si a.t g = 1. Data points: Monte Carlo result for A'^ = 10. Solid line: Best 
free variational approximation. 

FIG. 4. Eigenvalue distribution p(A) generated by S*! at (7 = 1 in Monte Carlo simulations for A'^ = 20 (solid line), A = 10 
(dotted), A = 3 (dashed). 

FIG. 5. Eigenvalue distribution p(A) generated by Si at (? = 4. Data points: Monte Carlo result for N — 10. Solid line: 
Best free variational approximation, namely in the variables A/i + M2, Mi — M2. Dotted: Free approximation in the original 
variables Mi, M2. 

FIG. 6. Eigenvalue distribution p(A) generated by Si'^''. Data points: Monte Carlo result for N — 10. Solid line: Best free 
variational approximation. 

FIG. 7. Eigenvalue distribution p(A) generated by S2 a.t g = 2/5. Data points: Monte Carlo result for N — 40. Solid line: 
Best free variational approximation. 

FIG. 8. Eigenvalue distribution p(A) generated by Sd2- Monte Carlo results for A = 10 are given in the cases D — 2 
(dash-dotted), D = 3 (dashed) and D = 10 (dotted). Solid line: Free variational approximation in the original variables Mi. 

FIG. 9. Eigenvalue distribution p(A) generated by Sd4- Monte Carlo results for A = 10 are given in the cases D — 2 
(dash-dotted), D — 3 (dashed) and D = 10 (dotted). Solid line: Free variational approximation in the original variables A/i. 
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FIG. 10. Eigenvalue distribution p(A) generated by Sqi at g = 0.1. Data points: Monte Carlo result for A'' = 10 and L = 10. 
Solid line: Free variational approximation in the Fourier basis. Dotted: Free variational approximation in the original basis. 
The former is the better approximation of the two. 

FIG. 11. Eigenvalue distribution p(A) generated by Sqi at g = 1. Data points: Monte Carlo result for AT = 10 and L = 10. 
Solid line: Free variational approximation in the Fourier basis. Dotted: Free variational approximation in the original basis. 
The former is the better approximation of the two. 

FIG. 12. Eigenvalue distribution p(A) generated by Sqi at g — 10. Data points: Monte Carlo result for A'" = 10 and L = 10. 
Solid line: Free variational approximation in the Fourier basis. Dotted: Free variational approximation in the original basis. 
The latter is the better approximation of the two. 

FIG. 13. Eigenvalue distribution p(A) generated by Sqi in the continuum limit, namely for g = 1 (for the definition of g, see 
text). Dotted: Exact result. Solid line: Free variational approximation in the Fourier basis. 

FIG. 14. Eigenvalue distribution p(A) generated by Sq2 at g = 0.1. Data points: Monte Carlo result for AT = 10 and L = 10. 
Solid line: Free variational approximation in the Fourier basis. 

FIG. 15. Eigonvahie distribution p(A) generated by Sq2 at g = 1. Data points: Monte Carlo result for N = W and L = 10. 
Solid line: Free variational approximation in the Fourier basis. 

FIG. 16. Eigenvalue distribution p(A) generated by Sq2 at g = 10. Data points: Monte Carlo result for N = 10 and L = 10. 
Solid line: Free variational approximation in the Fourier basis. 
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